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We calculate the mean-field thermodynamics of a spherically trapped Fermi gas with unequal 
spin populations in the unitarity limit, comparing results from the Bogoliubov-de Gennes equations 
and the local density approximation. We follow the usual mean-field decoupling in deriving the 
Bogoliubov-de Gennes equations and set up an efficient and accurate method for solving these equa- 
tions. In the local density approximation we consider locally homogeneous solutions, with a slowly 
varying order parameter. With a large particle number these two approximation schemes give rise 
to essentially the same results for various thermodynamic quantities, including the density profiles. 
This excellent agreement strongly indicates that the small oscillation of order parameters near the 
edge of trap, sometimes interpreted as spatially inhomogeneous Fulde-Ferrell-Larkin-Ovchinnikov 
states in previous studies of Bogoliubov-de Gennes equations, is a finite size effect. We find that 
a bimodal structure emerges in the density profile of the minority spin state at finite temperature, 
as observed in experiments. The superfiuid transition temperature as a function of the population 
imbalance is determined, and is shown to be consistent with recent experimental measurements. 
The temperature dependence of the equation of state is discussed. 

PACS numbers: 03.75.Hh, 03.75.Ss, 05.30.Fk 



I. INTRODUCTION 

There has been considerable recent experimental 
progress in creating strongly interacting ultra-cold 
atomic Fermi gases. A primary tool for the manipula- 
tion of these systems is the use of Feshbach resonances, 
through which the magnitude and sign of the inter-atomic 
interaction can be tuned arbitrarily by an external mag- 
netic field. For a two-component {i.e., spin 1/2) Fermi 
gas with equal spin populations, it has been expected 
for some time that the system will undergo a smooth 
crossover from Bardeen-Cooper-Schrieffer (BCS) super- 
fluidity to a Bose-Einstein condensate (BEC) of tightly 
bound pairs. This scenario has now been confirmed un- 
ambiguously by some recent measurements on both dy- 
namical and thermodynamical properties El, [3, 0, 0, S, S, 

Since the population in each s pin state can also be 
adjusted with high accuracy [l^, [ill, [l2l, subtle 
question of particular interest is the ground state of 
a spin-polarized Fermi gas with different particle num- 
bers in the spin up and down states. As conventional 
BCS pairing requires an equal number of atoms for 
each spin component, exotic forms of pairing are nec- 
essary in order to accommodate a finite spin popula- 
tion imbalance. There are several scenarios suggested 
in the weakly coupling BCS limit for a uniform gas, 
including the spatially modulated Fulde-Ferrell-Larkin- 
Ovchinnikov (FFLO) st ate fisll , the breached pairing [l3| 
or Sarma superfiuidity |l5l. lid. [T^ . and phase separa- 
tion [18]. In the strong coupling, BCS-BEC crossover 
regime, a variety of mean-field phase diagrams have been 
also proposed [H, M, HH, [13, M, [H, [11, [1^, [IS- How- 



ever, no clear consensus on the true ground state of spin- 
polarized fermionic su perfi uidity has been reached as yet 

0, M, M, M, M, MM, 35, 361 

Recent investigations [S, ilfl, Till . [l3| on atomic ^Li 
gases with tunable population imbalance open up intrigu- 
ing possibilities for solving this long-standing problem. 
These experimental observations have attracted intense 
theoretical interest Js^, M, M, \M, \M, [13, [13, \M, \M, 
\M, [13, [H, \M, \M, \5M- We note that there is no firm 
experimental evidence for the various non-standard su- 
perfiuid states mentioned earlier, which involve homoge- 
neous spin-polarized environments. Various interesting 
phenomena have been demonstrated experimentally, in 
optical traps of different shapes and sizes: 

(A) A shell structure is observed in the density profiles 

by Zwierlein et al pjoj, with a bimodal distribution 
at finite temperature, suggesting an interior core of 
a BCS superfiuid phase with an outer shell of the 
normal component. As a result, the thermal wing 
may provide a direct route to thermometry, in an 
environment where temperature is often difficult to 
calibrate reliably. 

(B) With increasing spin-polarization, the gas shows a 

quantum phase transition from the superfiuid to 
normal state. Close to the broad Feshbach res- 
onance of ^Li at B ~ 833 G, a critical polar- 
ization Pc = 0.70(3) has been determined at low 
temperatures. Here, the relative polarization is 
p = [N^ - iV|)/(7V| + iVx) where N„ is the num- 
ber of spin up or down atoms. The value of Pc 
decreases with increased temperature. 

(C) A similar shell structure is also identified by Par- 
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tridge et al [I4I, in an experimental trap configura- 
tion with a very large aspect ratio. However, this 
bimodal structure disappears below a threshold po- 
larization of P, ~ 0.10. 

To make quantitative contact with the current ex- 
perimental findings, it is crucial to take into account 
the trapping potential that is necessary to prevent the 
atoms from escaping. Therefore, the theoretical analysis 
is more complicated. The simplest way to incorporate 
the effect of the trap is to use a local density approx- 
imation (LDA), where the system is treated locally as 
being homogeneous, with spatial variation included via 
a local chemical potential that includes the trap poten- 
tial. Although this method has been extensively used to 
study the densi ty p rofiles of s p in- p olarized Fermi gases 
[13, S, m, 53, EJ, Si, Si, SI, SH, Sil, its validity has 
never been thoroughly examined. 

Alternatively, within the mean-field approximation, 
one may adopt the Bogoliubov-de Gennes (BdG) equa- 
tions, which include the full spatially var ying trap po- 
tential from the outset [13, SI, S3, [13, Isif ! It has 
been claimed that the solution of BdG equations in- 
cludes FFLO states with spatially varying order param- 
eter [H, S3, HO, IHH . However, the only evidence for this 
statement is the observation that in such solutions the 
sign and amplitude of the order parameter or gap func- 
tion exhibit a small oscillation near the edge of traps. 

In this paper, we perform a comparative study of the 
thermodynamic properties of a trapped, spin-polarized 
Fermi gas, by using both the mean-field BdG equations 
and the LDA. While neither method is exact, since pair- 
ing fluctuations beyond the mean-fleld approximation are 
neglected, this comparison can give insight into the con- 
sequences of the different types of approximation in com- 
mon use. As the most interesting experiments were in the 
BCS-BEC crossover regime, we focus on the on-resonance 
situation. In this regime - sometimes called the unitary 
limit [131, the s-wave scattering length diverges. 

The purpose of this work is three-fold. First of all, we 
have developed an efficient and accurate hybrid method 
for solving the mean-field BdG equations, based on the 
combined use of a mode expansion in a finite basis for 
low- lying states, together with a semiclassical approxi- 
mation for the highly excited modes beyond a suitably 
chosen energy cut-off. The cut-off is then varied to check 
the accuracy of the hybrid approach. As a consequence, 
we are able to consider a Fermi gas with a large number of 
atoms (~ 10^) that is of the same order as in experiment. 

The second purpose is to use this gain in efficiency to 
perform a detailed check on the accuracy of the LDA 
description in comparison with the BdG results. It is 
worth noting that, unHke previous theoretical work, we 
do not include the Hartree term in the BdG equations, 
since this should be unitarity limited in the BCS-BEC 
crossover regime [53|. For a sufficiently large number of 
atoms, we find an excellent agreement between these two 
approximation schemes. In particular, small oscillations 
of the order parameter at the edge of traps which were 



reported previously, tend to vanish as the number of par- 
ticles rises. Therefore, we interpret this as a finite size 
effect, rather than the appearance of a spatially modu- 
lated FFLO state. 

Finally, various thermodynamical quantities are calcu- 
lated. The observed bimodal structure in the density dis- 
tribution of the minority spin component is reproduced 
theoretically. The transition temperature is determined 
as a function of the spin population imbalance, and is 
found to qualitatively match the available experimental 
data. The temperature dependence of the equation of 
state is also discussed. 

The paper is organized as follows. In the next section, 
we present the theoretical model for a spin-polarized, 
trapped Fermi gas. In Sec. Ill and Sec. IV, we ex- 
plain the LDA formalism, and then describe in detail how 
we solve the BdG equations. The relationship between 
these two methods is explored. In Sec. V, a detailed 
comparison between LDA and BdG calculations is per- 
formed. Results for various thermodynamical quantities 
are shown. Their dependence on the population imbal- 
ance and on the temperature are studied. Sec. VI gives 
our conclusions and some final remarks. 



II. MODELS 

We consider a spherically trapped spin-polarized Fermi 
gas at the BCS-BEC crossover point, as found near a 
molecular Feshbach resonance. In general, a system like 
this requires a detailed consideration of molecule forma- 
tion channels [Hi, but near a broad Feshbach resonance, 
the bound molecular state has a very low population[55|. 
Accordingly, the system can be described approximately 
by a single channel HamiltonianjH^, US], 



2m 



(r) 



+ U J d^r*+ (r) (r) (r) (r) , (2.1) 

where the pseudospins a denote the two hyperfine 

states, and (r) is the Fermi field operator that anni- 
hilates an atom at position r in the spin a state. The 
number of total atoms is N ~ + Ni. Two different 
chemical potentials, ^1,1 — ^i±S^, are introduced to take 
into account the population imbalance SN = N-^ — N^, 
V (r) = muj'^r'^/2 is the isotropic harmonic trapping po- 
tential with the oscillation frequency u, and U is the bare 
inter-atomic interaction strength. We now describe the 
LDA and BdG theories. 



III. LOCAL DENSITY APPROXIMATION 

If the number of particles becomes very large, it is 
natural to assume that the gas can be divided into many 
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locally uniform sub-systems with a local chemical po- 
tential [l^l- Then, within the LDA, the trap terms in 
the Hamiltonian Eq. (|2.ip are absorbed into the chemi- 
cal potential, so that we have effective space-dependent 
chemical potentials. 



MT (r) = A*T - ^ (r) , 
A^i (r) = Hi-V{r). 



(3.1) 



Note that the local chemical potential difference 5^{v) = 
[^1 (r) — /i| (r)] /2 = 5(1 is always a constant, but the av- 
erage /i (r) = (r) -I- Hi (r)] /2 decreases parabolically 
away from the center of trap. 



A. Effective Hamiltonian 

If the global potentials and are fixed, we can 
consider a locally uniform Fermi gas in a cell at position 
r with local chemical potentials /i| (r) and (r), whose 
Hamiltonian takes the form, 



H{v) 



kcr 



2m 



- /"<T (r) 



kk'p 



kf Cp_j^j^Cp_k'J,Ck'| . 



(3.2) 



Here Cka represents the annihilation operator for an atom 
with kinetic energy fi?k'^ / (2m). For simplicity we re- 
strict ourself to homogeneous superfiuid states for the 
locally uniform cell, i.e., the local order parameter has 
zero center-of-mass momentum. By taking the mean- 
field approximation, an order parameter of Cooper pairs 
A (r) = X)k ('^klC-kf) is therefore introduced, whose 
value depends on position owing to the spatial depen- 
dence of (i-\ (r) and (r). The local Hamiltonian (|3.2p 
then becomes, 



Wm/ (r) 



E 

kcr 



2m 



- M<T (r) 



- A(r) 



2^ [CkiC-kT 



'-kcrCktr 

h.c] . 



A2 (r) 



U 



(3.3) 



We have neglected the Hartree terms U ( Cj^tCkt ) and 



[/ X^k \'^ki'^kiy in this mean-field factorization. Their 
absence is owing to the following reasons: 

• The use of contact interactions leads to an unphys- 
ical ultraviolet divergence, and requires a renor- 
maHzation that expresses the bare parameter U 
in terms of the observed or renormalized value 
[A-Kh^a/m) , i.e., 



1 

U 



?2ek' 

k 



(3.4) 



where a is the back-ground s-wave scattering length 
between atoms, and ek = fi^k"^ / (2m). Generically, 
this renormalization requires an infinitely small 
bare parameter, in order to compensate the ultravi- 
olet divergence in the summation l/2ek. There- 
fore, strictly speaking, within a mean-field approxi- 
mation the Hartree terms should vanish identically. 

• For weak couplings, one may indeed obtain Hartree 
terms like (47rft^a/m)nx_x. With renormalization, 
these corrections are beyond mean-field, and are 
effective only in the deep BCS limit. Towards the 
unitarity limit with increasing scattering length, 
they are no longer the leading corrections and 
become even divergent. Higher order terms are 
needed in order to remove the divergence at uni- 
tarity. For example, one may use Fade approxima- 
tions in the equation of state^^. Thus, throughout 
the BCS-BEC crossover region, the neglect of the 
Hartree terms is not an unreasonable approxima- 
tion. 

The above mean-field Hamiltonian can be solved by the 
standard Bogoliubov transformation jHl|. The resulting 
mean-field thermodynamic potential has the form. 



^^m/ (r) 



(r) + y: 



A2 (r) 



2ek 



i5][ln/(-ii;k+)+ln/(-i?k-)], (3.5) 



f3 



where / (x) = [exp (/3x) + 1] ^ is the Fermi distribu- 
tion function (/3 = l/fc^T), and i?k± — -E-k ± SiJ,(r) 

are the quasiparticle energies with i?k — [Ck + (i")] 
and ^k = h?k'^/2m — /i(r). Given local potentials /i (r) 
and 6(1 (r) , we determine the value of order parameters 
A (r) by minimizing the thermodynamic potential, i.e., 
dQmf (r) /dA (r) = 0, or expHcitly, 



E 

k 



1 

2^ 



l-/(i?k+)-/(i?k-) 



(3.6) 



We note that a non-BCS superfiuid solution, the so-called 
Sarma state [H, [H, may arise in solving the gap 
equation p.6p . However, on the BCS side such solu- 
tion suffers from the instabilities with respect to either 
the phase separation or a finite-momentum paired FFLO 
phase Further, the Sarma state is not energetically 
favorable [131, and thereby will be discarded automati- 
cally in the numerical calculations. 



B. Thermodynamic Quantities 

Once the local order parameter is fixed, it is straight- 
forward to calculate the various thermodynamic quan- 
tities. The local particle densities are calculated ac- 
cording to, n| (r) = —dflmf (r) (r) and ni (r) = 
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-dn^nf (r) /dfii (r), or, 

k 



2Ek 



2Eu 



fiE^-) + ^^fi-E^+) 

/(i?k+) + ^^^/(-£;k-) 



2Ei, 



(3.7) 



while the entropy and the energy are determined, respec- 
tively, by 



S{r) = -ks E [/(Ska)ln/(£;kc 

k,Q=± 

+ /(-Ska)ln/,(-£;ka)], 



(3.8) 



and. 



S (r) = n^f (r) + T5(r) + n^n^ (r) + ^^n^ (r) . (3.9) 
Then, the integration over the whole space gives rise to, 

N (/i, 6fi) = J d^r [ni (r) + (r)] , 

SNi^l,6^l) = I rf^rK {v)-n^ (r)] , (3.10) 



and S = Jd^rS{r), E = ^d^rE{Y). The global chem- 
ical potentials ^ and S^jl should be adjusted to satisfy 
N{fi,Sfj.) = N, and SN{pL,6fj,) = 5N. The numerical 
calculation thereby involves an iterative procedure. 

We note that, on physical grounds, a general pic- 
ture can be drawn for the density profiles Near 
the center of the trap, where is small compared to 
fi (r) , the densities of the two spin states are forced 
to be equal and we have a BCS core extended up to 
a radius Rbcs- Outside this radius a normal state 
is more favorable than a superfluid phase. At zero 
temperature, the Thomas-Fermi radius of the minority 
(spin down atoms) and majority (spin up atoms) are 

given by R^^p = [2 (fi — S^i) / (^muj'^)]^^'^ and R^^p — 

1 /2 

[2 (p, + Sfi) I (mw^)] , respectively, as we neglect the 
interactions between the two components in the normal 
state. 

We finally remark that this entire approach is less 
accurate at Feshbach resonance, and especially on the 
BEC side of resonance, where it becomes essential to 
include quantum fluctuations beyond the mean-fleld 
approximationfg^l ■ 



IV. THE BOGOLIUBOV-DE GENNES 
MEAN-FIELD THEORY 

We next consider the Bogoliubov-de Gennes theory of 
an inhomogeneous Fermi gas, starting from the Heisen- 
berg equation of motion of Hamiltonian l|2.ip for (r, t) 
and *x (r,t): 



ih, 
ih 



9* 



dt 
dt 



2m 



2m 



(4.1) 



A. Mean-field approximation 

Within the mean-field approximation as before, we 
replace the terms J7*^*x*| and U^^^i^^ by their 

respective mean-field approximations J7*x^|4'^ = 
-A(r)*+ + C/nx(r)^'| and C/^'x*t*:[ = -A(r)^'+ + 
f7rt|(r)\l/x, where we define a gap function A(r) = 
-t7(*X*l) and n<x(r) = The Hartree terms 

C/rtcr(r) are infinitely small in the mean-field treatment 
due to the regularization of the bare interaction ?7 ^ 0, 
as discussed in greater detail below. We keep them in 
the derivation at the moment for clarity. The above de- 
coupling thus leads to. 



dt 

7 ± 

dt 



[Hf-Mi] *i + A(r)*+, 



(4.2) 



where H% = -fi2y2/ (2m) + V{r) + Uug. (r). We solve 
the equation of motion by the insertion of the standard 
BogoHubov transformation: 



-iEjyt/h/ _|_ ^+ g-iBj^t/n 



(4.3) 

This yields the well-known BdG equations for the Bo- 
goliubov quasiparticle wave functions uja- (r) and Vja- (r) 
with excitation energies Ej^, 



K - A(r) 
A*(r) + 





Uja (r) 


- Eja- 


■ Uj„ (r) " 




. Vja (r) _ 


. Vja (r) _ 



(4.4) 



where Uja- (r) and Vja (r) are normalized by 
/ o?r( \ujc (r)l^ + \vja (r)l^) = 1. The number den- 
sities of different hyperfine states (r) = (^^'|^^'x^ and 

n| (r) = Z^'+^'x V and the BCS Cooper-pair condensate 
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A(r) = — can then be written as, 

(4.5) 

where the statistical averages {Cj^Cja) — f (Eja) and 
{cjaC^^) — f {—Eja) have been used. 

The solutions of the BdG equations contain both posi- 
tive and negative excitation energies. Thus, to avoid dou- 
ble counting, a factor of half appears in the summation in 
Eq. I|4.5p . Furthermore, the presence of the chemical po- 
tential difference breaks the particle-hole symmetry and 
therefore leads to different quasiparticle wave functions 
for the two components. One can easily identify that 
there is a one to one correspondence between the solu- 
tion for the spin up and spin down energy levels, i.e., 



and 



Eja ^ ^Ejg-, 



and 



Uja (r) 
Vja (r) 



(r) 



(4.6) 



(4.7) 





A(r) ■ 




■ Uj (r) ■ 




Uj (r) 


A*(r) 


-HI + m 








. (r) _ 



By exploiting this symmetry of the BdG equations, there- 
fore, we need to solve the BdG equations for the spin 
up part only. This has the following form after remov- 
ing the spin index, i.e., we let Uj (r) — Uj-\ (r) and 
«i (r) = fjT to give: 



(4., 



B. Hybrid BdG Technique 

We now wish to address the issue of how to use the 
BdG equations in a practical numerical application. Ac- 
cordingly, the density distributions and the gap function 
can be rewritten as, 

n^{v) = EK(r)|V(i?,), 

j 

3 

A(r) = C/Eu,(r)z;*(r) /(£;,). (4.9) 



Eqs. I|4.8p and l|4.9p can then be solved self-consistently, 
with the constraints that 

N 5/i) = / d^r [nt (r) -f (r)] =iV, (4.10) 



6N {fi, Sfi) = / dh [n^ (r) - (r)] ^SN. (4.11) 



where. 



In any practical calculation, due to computational lim- 
itations, one has to truncate the summation over the 
quasiparticle energy levels. For this purpose, we intro- 
duce a hybrid strategy by introducing a high-energy cut- 
off Ec, above which the local density approximation may 
be adopted[60| for sufficiently high-lying states. Follow- 
ing this approach, we then have, 

n|(r) = n|,d(r)+n|^c(r), 
= "i,d(i")+nix(r), 
A(r) = A<i(r) + A,(r), (4.12) 

\Ej\<E, 

r.T.c(r) = J2 K(r)lV(i?,), 

\Ej\>Ea 

n^,,(r) = Yl \vAr)\'fi-Ej), (4.13) 

\Ej\>Ea 

A,(r) ^ U Y %(r)«*(r)/(i?,), 

\Ej\<E, 

A,(r) = U J2 u,{v)v*{v)f{E,). (4.14) 

\E,\>E^ 

We consider below separately the contributions from the 
quasi-continuous high-lying states {\Ej\ > Ec) and dis- 
crete low-energy states {\Ej\ < Ec). This allows us to 
take into account the spatial variation of the low-lying 
trapped quasiparticle wave functions, without having to 
treat all high-energy states in this formalism. 

C. LDA for high-lying states 



and 



For the BdG equations l|4.8p , the local density approx- 
imation is the leading order of a semiclassical approxi- 
mation and amounts to setting 

Uj (r) u (k, r) exp [ik • r] , 
Vj (r) — > V (k, r) exp [ik • r] , 

Ej -> £;(k), (4.15) 

where u (k, r) and v (k, r) are normalized by \u (k, r)|^ -|- 
\v (k, r)|" = 1, and the level index "j" has now been re- 
placed by a wave vector k. So the equations in l|4.8p are 
reduced to the algebraic form 
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(k,r) 



A(r) 



A*(r) 



m (k,r)+/i4 



£;(k) 



(4.16) 

where the quasi-classical single particle Hamiltonian is 

h^k^l {2m 



K (k,r) 



V (r) + Un^ (r) 



(4.17) 



We obtain two branches of the excitation spectra, 
£;(k,+) = Ei,-Sfi-U[n^{r)-ni{r)]/2 , andE{k,-) = 
Ei^ + Sfj. + [/[n-f (r) — ni{r)]/2, which may be interpreted 
as the particle and hole contributions, respectively. Here 

Ej, = [^l + A2(r)] and ^k = h^ky2m + F(r) - ^l + 
Un{v)/2. Note that consistent with the definition in Sec. 
Ill, we have reversed the sign of the excitation spec- 
trum of the hole branch, that is, for the particle branch 
£;(k) = +i;(k,+), while for holes E{\<i) = -E{k,-). 
The eigenfunctions of the two branch solutions are, re- 
spectively. 



1 



< = 2 I 1 



Sk 

in 

2 Sk 
A(r) 
2E^ ' 



(4.18) 



and 



2 1 ?k 

ut — - 1 

2 V Sk 



UkWk 



A(r) 

2£;k ^ 



(4.19) 



Thus, above the energy cut-off, the quasi-continuous con- 
tribution of high-lying states to the density profiles and 
the gap function can be obtained by, 



£;(k,+)>£;c 



^ fi-Ejk,-)] 

^ 2 V £^k 



n\,c (r) 



E 

£;(k,+)>£;c 



/h^(k,+)] A _ ik 
2 V Sk 



£;(k,-)>£;c 



and 



Ac(r) = E 



£;(k,+)>_Bc 



A(r) 
2i;k 



/[i?(k,+)] 



E (4.21) 

£;(k,-)>£;c 



It is worth noting that if we reduce the energy cut-off 
Ec to zero, we recover the LDA expressions for the den- 
sity profiles in Sec. Ill (see, for example, Eq. (|3.7p ). 
Moreover, Eq. (|4.2ip reduces to the LDA gap equation 
(011). 

At the other extreme, for a sufficiently large energy 
cut-off {(iEc » 1), we may discard the Fermi distribution 
function in Eqs. I|4.20p - I|4.2ip . As a result we have the 
following simplified gap equations for the above cut-off 
LDA contributions: 



and 



"T,c(r) 
"ix(r) 

Ac(r) 



in 

E^ 



E - 

^ 2 

E(k,-)>£;c 

E -fi-- 



£;(k,-|-)>_Ec 



^ E 

£;(k 



A(r) 



2E^ 



(4.22) 



(4.23) 



D. BdG equations for low-lying states 

Let us now turn to the low-lying states by solving the 
BdG equations l|4.8p . As we consider a spherical trap, it 
is convenient to label the Bogoliubov quasiparticle wave 
functions Uj (r) and Vj (r) in terms of the usual quantum 
numbers j — {nZm}, and write. 



u, r 



Unl (r) 



(4.24) 



Here u„; (r) /r and Vni (r) jr are the standard radial wave 
functions, and Yim {d, (f) is the spherical harmonic func- 
tion. The BdG equations are then given by. 



A(r) 



A(r) 





Unl 


— E„i 


Unl 




Vnl 


Vnl 






(4 



where 

K (0 



2m 



1(1 + 1) 



+ V{r) + Una{r) (4.26) 



is the single particle Hamiltonian in the I sector. 
We solve these equations by expanding Uni (r) and 
Vnl with respect to the eigenfunctions 4'ai (r) of a 
3D harmonic oscillator radial Hamiltonian, Hose (0 = 
-nV {2m) [d^/dr^ + I {I + 1) /r^] + V {r) . These have 
energy eigenvalues eai — {2a + I + 3/2)fku, and the re- 
sulting expansion is: 



Unl {r) = J2^ni(l^°^i (r) 
Vnl (r) = ^B'^i(j)c,i (r) 



(4.27) 
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The problem is then converted to obtain the eigenvalues 
and eigenstate of a symmetric matrix, 



Aa/3 



T 



i 



E„ 



4" 



where, we have defined Caia = — A'cr, and: 
Aa/3 = / dr(j)ai (r) A (r) (r) , 



(4.28) 



M. 



Q/3 



(r) Uni {r)4>i3i (r) 



(4.29) 



We note that the renormalization condition Uni {r) and 
Vni (r), / dr [w^; (r) + (r)] = 1, is strictly satisfied, 

since Y.a i^nif + i^nif = 1- Once Uni (r) and u„; (r) 
are obtained, we calculate the gap equation for the low- 
lying states {\Eni\ < E^) and the corresponding number 
equations: 



A4r) 

n^,d (r) 



2Z + 1 

^^^2 M«i (r) w„i (r) / (£■„;) , 



E 

nl 

E 



2? 



475-7-2 
21 + 1 

47^7-2 



vli {r)f{~E^i). 



(4.30) 



E. Regularization of the bare interaction U 

We now must replace the bare interaction U by the 
corresponding s-wave scattering length, using standard 
renormalization techniques. The combination of expres- 
sions (|4.23p and l|4.3Qp gives the full gap equation, 

\Ej\<E^ E{k.-)>Ec 

(4.31) 

which is formally ultraviolet divergent due to the use of 
the contact potential. However, the form of the sec- 
ond term on the right side of the equation suggests a 
simple regularization procedure. We substitute 1/U — 

{Anfi^a/rri) — X^k^/^^k into the above equation and 
obtain, 

;A(r) = n,{v)v*{v)f{E,) 

\E,\<E^ 



A(r) 



A(r) 



2ew E/ 2E\^ 

k £;(k,-)>_Bc 



(4.32) 



Thus we may rewrite the gap equation in terms of an 
effective coupling constant Ueff (r), i.e., 



A(r) = U,ff (r) 



ff ^ ' ^ 

\E,\<E, 



u,ir)v*{v)fiE,), (4.33) 



where we have introduced an effective coupling constant 
Ueff (r) defined so that: 



^ 2ek 
k " 



E — 



E{k,-)>Ec 

(4.34) 

The ultraviolet divergence now cancels in the bracketed 
expression. The effective coupling constant Ueff (r) will 
depend on the cut-off energy. However, the resulting gap 
in Eq. I|4.33p is essentially cut-off independent. 

The use of the uniform regularization relation 1/U = 

(^Anh^a/nij — X^k^Z-^^k leads to an infinitely small 
bare interaction coupling. One therefore has to replace 
U by zero anywhere if there is no ultraviolet divergence 
in the summations. As mentioned earlier, this replace- 
ment is the proper treatment within mean-field theory. 
Certainly, this procedure neglects the Hartree correction, 
which is of importance in the deep BCS regime. How- 
ever, around the unitarity regime of interest here, the 
usual expression for the Hartree correction becomes di- 
vergent, and requires a more rigorous theoretical treat- 
ment which shows that it is no longer significant (58|. 
Consistent with this treatment, we note that these mean- 
field Hartree shifts are not observed experimentall y in the 
energy spectra in the BCS-BEC crossover regime [53|- In 
other words, the Hartree terms should be unitarity lim- 
ited at crossover. 

It is important to point that in principle, the regu- 
larization procedure proposed above is equivalent to the 
use of a pseudopotential, as suggested by Bruun and co- 
workers [61|. However, the pseudopotential regulariza- 
tion involves a calculation of the regular part of the Green 
function associated with the single particle Hamiltonian 
Hs and is numerically inefficient. Alternative simplified 
regularization procedures have also been introduced by 
Bulgac and Yu [62], and Grasso and Urban Our 
prescription l|4.34p may be regarded as a formal improve- 
ment of these regularization procedures. 



F. Summary of BdG formalism 

We now summarize the BdG formalism, by converting 
the summation over the momentum k in the high-lying 
levels to a continuous integral of the energy. 

We find that the total spin densities are given by: 

-T(r) = E 

|£;„il<£;c 

+ / den|,c (e,i") , 



, , 2/ + 1 2 

-i (r) - E 1^ ' 

\E„i\<E, 



vli {r)n-Er.i) 



+ / c^e'^i.c (e, r) , 



(4.35) 
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with a modified gap equation of: 



A(r) 
Ueff (r) 



^^^Unl (r) Vnl (r) f (Enl) . 



\E„i\<Ea 



(4.36) 



The above-cutoff contributions are given by: 



riTx(e,r) 



ni,c (r) 



e + 6^ 



1/2 



e — S^j, 



V(e-<5M)-A^(r) 



1/2 



(4.37) 



and the value under the square root is understood to be 
non-negative. Moreover, in an integral form the effective 
coupHng takes the form, 



1 



Ueff{r) Anfi^a 2t:^ Air^h? J e, 



def (e, r) , 



where. 



(4.38) 



1/2 



/(e,r) 



1/2 



v/(e-5/i)-A2(r) 



V(e-<5M)-A2(r)+Ai-y 



- 1 



(4.39) 



The radial wavefunctions in Eqs. (|4.35p and l|4.36p are 
calculated by solving the eigenvalue problem l|4.28p . As 
the matrix involves the gap function, a self-consistent 
iterative procedure is necessary. For a given number of 
atoms {N = N-^ + Ni and SN — N-^ — ^j,), temperature 
and s-wave scattering length, we: 

1. start with the LDA results or a previously deter- 
mined better estimate for A (r), 

2. solve Eq. I|4.38p for the effective coupHng constant. 



3. then solve Eq. (|4.28p for all the radial states up to 
the chosen energy cut-off to find m„; (r) and Vni (r), 
and 

4. finally determine an improved value for the gap 
function from Eq. I|4.36p . 



During the iteration, the density profiles (r) 
and ni (r) are updated. The chemical poten- 
tials and Sfi are also adjusted sHghtly in each 
iterative step to enforce the number-conservation 
condition that drAirr"^ [n^ (r) -I- ni (r)] —N and 
drAnr^ [n| (r) — ni (r)] —SN, when final conver- 
gence is reached. To make contact with the experimental 
observed density profiles (lol . [l^ , we calculate the axial 
and radial column densities, 



J —QC ^ 

niip) = f dzn^(Vp^T^), (4.40) 



and 



ni (z) = J 2Trpdpni (^\/ p^ + , 

nx (z) = ^ 2Tipdpn^ (^fp^Tz^) ■ (4.41) 



G. Entropy and energy 

Apart from the density profiles and gap function, we 
can also determine the entropy and total energy of the 
imbalanced Fermi gas, by using the expressions: 



-fcs ^ (2Z 1) [/ {Enl)\nf{Enl) 



-f{-Eni)\nf{-E„ 



and 



E^ j d^v\^[^+ {r)n%^a (r)] 



The energy can further be written as, 

m 



E = 



PlNi + PiNi - 



|A(r) 

U 



(4.42) 



(4.43) 



j 

fd^rY^ 

'J 1, 



2ek 



(4.44) 



where we have replaced the bare interaction U by the 
s-wave scattering length. The contribution of the high- 
energy part to the entropy is essentially zero. 

For the total energy, we must take into account both 
the low-lying states and the high-lying states. Therefore, 
we divide the energy into two parts E = Ed + Ec, where 



E, = 



(2^ + 1)^" 

|s„,t<£c 



drA 



TTT^ lAI 



fiE^i)- I drviiir) 





(4.45) 
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and 



= / d'V 



E 1 



E(k,+)>Ec 



£;(k,-)>£;c 



in 



i?(k,+) 



fk 



(4.46) 



By converting the summation into an integral, we obtain 
Ec — drAnr'^Ec (r), where 



Eciv) 



|A(r)r [fc,,i fce.2 

2 1 27r2 27r2 



A [£;,a (e, r) + ^c,2 (e, r)] | , (4.47) 



fcc.i = [V{E, + SiJ.)-A''ir)+^i-Vir) 
kc.,2 - [v/(-Ec-<5m)- A2(r)+/i-F(r) 



1/2 
1/2 



(4.48) 



and 



^^c,i(e,r) 



Ec,2 (e, r) 



y2m3/2 [£+ + ^ - y (r)] 



1/2 



E. 



e (e + (5/x) /2 

e + 6^i + E+ E++fi-V (r) 

\/2m3/2 +/^- V(r)]'/' 
ft3 A2(i.) 

(e - 5/i) /2 



e-Sfi + E^ E^+ fi-V (r) 



(4.49) 



(i.e., iV| = A''| = N/2). At zero temperature, a sim- 
ple LDA treatment of the harmonic potential leads to 

1 /3 

a Fermi energy and a Fermi temper- 

1/3 

ature Tp = {3N) hw/kB- Accordingly, the Thomas- 
Fermi (TF) radius of the gas is ttf — (24iV)^^^ a?io, 
and the central density for a single species is nTF(O) = 

In a uniform situation, the universality argument gives 
/i = ^ep in the unitary limit (H3|, where ep = h'^kp/{2m) 

1 /3 

with a Fermi wavelength kp = (STr^n) . The mean- 
field BCS theory predicts £, ~ 0.59. Therefore, a bet- 
ter mean-field estimate of the chemical potential and 
the TF radius for a trapped unitarity Fermi gas will be 

fJ-TF,unitarity — S^^'^Ep and TTP^unitarity = C"^^**?'TF- 

For most calculations, we use a number of total atoms 
N — 20000, which is one or two orders of magnitude 
smaller than in current experiments. We take a cut-off 
energy Ec = 64fiw, which is already large enough because 
of the high efficiency of our hybrid strategy. Typically 
we solve the BdG equations l|4.28p within the subspace 
n < JTrnax = 72 and I < /max — 120. The value of rimax 
and /max is determined in such a way that the subspace 
contains all the energy levels below Ec- The calculations 
usually take a few hours for a single run in a single-core 
computer with 3.0GHz CPU, for a given set of parameters 
T, N, SN, and iTrh^a/m. Further increase of the number 
of total atoms to 10^ or 10^ is possible, but very time- 
consuming. 

Below we will present some numerical results. In par- 
ticular we will examine the validity of the LDA at zero 
temperature and at finite temperature. We analyse some 
previous suggestions of the apparent appearance of the 
FFLO states in the mean-field BdG theory. Finally, 
we will investigate the thermodynamical behavior of the 
spin-polarized Fermi gas. 



LDA versus BdG 



where E± = y/{e ± Sfi) - A^{r). 



V. NUMERICAL RESULTS AND DISCUSSIONS 

To be concrete, we will focus on the on-resonance (uni- 
tarity) situation in our numerical calculation, in which 
the s-wave scattering length goes to infinity. It will be 
convenient to use "trap units", i.e., 



m = uj ^ h = ks = 1. 



(5.1) 



Therefore, the length and energy will be measured in 

1/2 

units of the harmonic oscillator length a^o — [fi/ {mu))] 
and huj, respectively. The temperature is then taken in 
units oihuj/kB. It is also illustrative to define some char- 
acteristic scales, considering a spherically trapped ideal 
Fermi gas with equal populations in two hyperfine states 



We present in Fig. 1 the density profiles of 
two spin states at temperatures T = O.OlTp and 
T — 0.20Tp for different population imbalances P = 
{N-^ - Ni) / {N^ + Ni) = 0.23, 0.48 and 0.84 as indi- 
cated. The results from the BdG and LDA approaches 
are plotted using solid lines and dashed lines, respec- 
tively. There is an apparent phase separation phe- 
nomenon, with a superfiuid inner core and normal shell 
outside, which is in consistent with the recent experimen- 
tal observation by Zwierlein et al. |il0|] and Partridge et 
al. [i3|. Particularly, for P = 0.84 at T = 0.20rF the 
minority (spin down) profile is enhanced at center, which 
in turn induces a slight decrease of the central density of 
the majority component. The appearance of a dense cen- 
tral feature in the minority spin profile agrees well with 
the on-resonance measurement reported by Zwierlein et 
al. [l3|. It clearly resembles the bimodal structure in the 
density distribution of a BEG. 
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Figure 1: (color online). Density profiles of spin up and spin down atoms for a spin-polarized unitary Fermi gas at different 
population imbalances and temperatures, as indicated in the figures. The number of total atoms is A'' = 20000. The density 
profiles are normalized by the Thomas-Fermi center density of an ideal symmetric Fermi gas with the same number of total 
atoms nTF{0) — {24N)^^'^ aj^^ / (Stt^), while the length is renormalized by the corresponding Thomas-Fermi radius ttf = 
(24A'')^'"' aho- The solid lines and dashed lines refer to the BdG results and LDA results, respectively. 



For all the spin polarizations considered, we find rea- 
sonable agreement between these two methods at the cho- 
sen total number of atoms N — 20000. As we shall see 
below, the agreement persists in various thermodynam- 
ical quantities, such as the chemical potential, entropy 
and total energy. >From the density profiles, the agree- 
ment is excellent at a small or intermediate population 
imbalance. The difference between the BdG and LDA 
prediction tends to be smaller as the temperature in- 
creases. For a large population imbalance, however, the 
agreement becomes worse. This can be understood from 
the corresponding gap functions as given in Fig. 2. The 
gap function in LDA experiences a sudden decrease at 
the superfiuid-normal interface. With increasing popula- 
tion imbalance, the drop is much more apparent, and ac- 
cordingly, the BdG gap function show a very pronounced 
oscillation behavior. Therefore, the deviation of the BdG 
density profiles from the LDA predictions becomes larger. 

It is important to note that only the axial column den- 
sity or radial column density can be measured by the ab- 
sorption imaging technique in the experiment. In Figs. 
3a and 3b, we plot respectively the axial column profile 
and radial column profile at T = O.OlTp for the imbal- 
ance P = 0.84. The minor difference between BdG and 
LDA shown in the three-dimensional density profiles is 
essentially washed out. 

B. FFLO state at the superfluid-normal interface? 

The consistency between BdG and LDA treatments 
reported here is in sharp contrast with some previous 



1.0 




0.0 0.2 0.4 0.6 0.8 1.0 

TF 

Figure 2: (color online). Gap functions at different temper- 
atures T = O.OITf (a) and T = 0.20Tf (b) for various pop- 
ulation imbalances as labeled. The number of total atoms is 
A'^ = 20000. The solid lines are the BdG predictions, while 
the dashed lines are the LDA results. The value of gap is 
renormalized by the non-interacting Fermi energy of an ideal 
symmetric Fermi gas Ef = (SA'^)^'^^ huj. Note that the small 
oscillation in the gap function at the superfiuid-normal inter- 
face at low temperature T = O.OITf vanish when the temper- 
ature becomes high enough. 
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Figure 3: (color online). Radial (a) and axial (b) column 
density profiles for iV = 20000 and P = 0.84 at T = O.OlTp. 
The difference between the BdG results (solid lines) and LDA 
results (dashed lines) becomes extremely small due to the 
integration over the extra dimensions. 
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Figure 4: (color online). Dependence of the gap function on 
the number of total atoms at T = O.OITf and P = 0.48. The 
solid line is the LDA result, and the others are BdG results 
with N = lO'^ (dash-dotted line), iV = 20000 (dotted line), 
and = 4000 (dashed line). 



studies |48l . l49l . ISOl . l5ll |. where a notable discrepancy of 
BdG and LDA results is found. In those studies the small 
oscillation of the gap functions at the superfiuid-normal 
interface is interpreted as the appearance of a spatially 
modulated FFLO state. Thus, the discrepancy is ex- 
plained as due to the breakdown of LDA approximation 
for FFLO states. From our results, the oscillation in the 
order parameters at the interface appears to be a finite 



size effect. This idea is supported by the observation that 
the BdG formalism naturally reduces to that of LDA if we 
set the cut-off energy Ec to zero as we mentioned earlier. 
On the other hand, as shown in Fig. 4, if we increase the 
number of total particles, the oscillation behavior of gap 
functions becomes gradually weaker. We thus infer that 
the oscillation will vanish finally in the limit of sufficient 
large number of atoms. 

To understand the discrepancy of BdG and LDA ap- 
proaches found in previous studies [isl . [49l . Isol . HH, sev- 
eral remarks may be in order. First, in these studies, 
the mean-field Hartree terms, i.e., Uni (r) and Un^ (r), 
appear in the decoupHng of the interaction Hamilto- 
nian in the BdG theory. However, the Hartree terms 
is absent in the corresponding LDA treatment. These 
terms cannot survive in the regularization procedure 
of the bare interaction U but are incorrectly included 
in Refs. Hi, Si, 113, [H as (47r ft^a/m) (r) and 
(Airfi^a/rnj (r), respectively. 

We note that the absence of the mean-field shift due 
to Hartree terms in the strongly interacting BCS-BEC 
crossover regime is already una mbig uously demonstrated 
experimental by Gupta et al. [53|- We conclude that 
there are three main reasons for the differences between 
the conclusion we find here that the two approaches are 
compatible, as opposed to earlier conclusions to the con- 
trary: 

1. The incorrect inclusion of Hartree terms in one ap- 
proach but not in the other, is the most Hkely rea- 
son for the discrepancy of BdG and LDA results 
shown in these previous works. 

2. The accuracy of numerical results depends crucially 
on the regularization procedure used to treat the 
bare interaction U . Proper treatment of regular- 
ization is essential. 

3. For a large population imbalance, the BdG equa- 
tions converge very slowly. Thus, a rigorous crite- 
rion is required to ensure complete convergence. 



C. Thermodynamic behavior of the imbalanced 
Fermi gas at unitarity 

We finally discuss the thermodynamics of the gas at 
unitarity. In Figs. 5, 6, and 7, we graph the chemi- 
cal potential, total energy and entropy of the g 
function of the population imbalance at various temper- 
atures. Again, we find good agreement between the BdG 
results and LDA predictions. With increasing population 
imbalance, the chemical potential decreases. For a given 
temperature there is a critical imbalance, beyond which 
the Fermi gas transforms into a fully normal state. The 
decrease of chemical potential becomes very significant 
when the phase transition occurs. In contrast, as pop- 
ulation imbalance increases, the total energy increases. 
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Figure 5: (color online). Chemical potentials as a function of 
the population imbalance at two temperatures T = O.OITf 
and T = 0.20rF. 



t4J 




Figure 6: (color online). Total energy per particle as a 
function of the population imbalance P at two temperatures 
T = O.OITf and T = 0.20Tf. Inset shows the first or- 
der derivation with respect to the population imbalance at 
T = 0.20rF. The jump at P ~ 0.9 marks the phase transi- 
tion to the normal state. 



Figure 7: (color online). Entropy per particle as a function of 
the population imbalance at various temperatures. 
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Figure 8: (color online). Average order parameter defined in 
Eq. I|5.2p as a function of the population imbalance at three 
temperatures T = O.OITf, T = 0.20rF, and T = 0.25rF. 
The position where the average order parameter vanishes de- 
termines the critical population imbalance towards the normal 
state at a fixed temperature. 



The impact of the phase transition on the total energy 
can scarcely be identified since the transition is smooth. 

As shown in the inset of Fig. 6, it can be exhibited 
clearly in the first order derivative of the energy with 
respect to the population imbalance, resembling the be- 
havior of the specific heat as a function of temperature. 
The entropy, on the other hand, shows a non-monotonic 
dependence at a temperature 0.20Tp. We identify this 
peak position as the phase transition point. 

To determine accurately the critical population imbal- 
ance as a function of temperature or vice versa, we plot 
in Fig. 8 an averaged order parameter at different tem- 
perature, which is defined via. 



N 



1/2 



(5.2) 



The condition ^ave = therefore determines the critical 
population imbalance Pc at a fixed temperature. The re- 
sulting Pc from the BdG calculation depends, of course. 



on the number of total atoms. To remove this depen- 
dence, we present the LDA prediction for critical tem- 
perature in Fig. 9a. To compare with the experiments, 
we also show four known experimental points in the fig- 
ure. The discrepancy between our theoretical predictions 
and the experimental data is most likely attributed to the 
strong pair fiuctuations beyond the mean-field, i.e., for a 
symmetric gas, the fiuctuation shifts Tc^bcs from 0.37Tf 
to around 0.27Tf [7], while at zero temperature, it reduce 
Pc,BCS from 1.0 to about 0.70 

Pair fluctuations must be taken into account in the 
strongly interacting unitarity regime for quantitative 
comparisons, as we have shown in earlier calculations 
[63 | without spin-polarization. This is certainly the 
most challenging problem in BCS-BEC crossover physics. 
Naively, we may linearly rescale the BCS critical tem- 
perature and population imbalance in such a way that 
both the theoretical Tc at P = and Pc at T = flts 
the experimental data. The third and fourth experimen- 
tal points of critical temperature at P = 0.56(3) and 
P = 0.59(3) now agree approximately with the re-scaled 
BCS Tc curve. 
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Figure 9: (color online), (a) BCS mean-field critical tem- 
perature as a function the population imbalance, determined 
from the LDA calculations (solid line). Open symbols are 
available experimental points from Zwierlein et al. [Tol . IT]] ] 
and Kinast et al. The dashed line shows a re-scaled BCS 
critical temperature as described in the text, (b) The critical 
entropy (the entropy at the critical temperature) against the 
population imbalance. 



For strongly interacting Fermi gases at low temper- 
ature, there is generally no reliable thermometry tech- 
nique. The bimodal structure in the density profile of the 
minority component rnay provide a useful method to de- 
termine temperature [10|. However its accuracy requires 
further theoretical investigations due to the strongly cor- 
related nature of the gas. Entropy is an alternative quan- 
tity that may be used to characterize the temperature in 
adiabatic passage experiments fg^l- In Fig. 9b, we show 
also the critical entropy of a trapped unitary gas against 
the population imbalance. The calculated dependence is 
essentially similar to that for temperature. 



tide number (~ 10^), the LDA appears to work very 
well. The discrepancy of BdG and LDA results reported 
in previous studies is attributed to the incorrect inclusion 
of a mean-field Hartree term in the BdG equations. We 
note, however, that the trap used in the current exper- 
iments is elongated in the axial direction. The spheri- 
cal trap considered in this work may be regarded as ap- 
proximately representative of the experimental setup by 
Zwierlein et al. [13, [ill with an aspect ratio approx- 
imately 5. The trap in the experiment by Partridge et 
al. [l3| is extremely anisotropic, and therefore the LDA 
description could break down. The solution of the BdG 
equations for such elongated systems is numerically in- 
tensive, and requires further investigation. 

Our derivation of the BdG formalism and the numeri- 
cal results with varying particle number in Fig. 4 strongly 
suggest that the calculated small oscillation in the order 
parameter at the superfiuid-normal interface arises from 
finite size effects. This is in marked contrast with the 
previous interpretation of this effect as due to a finite- 
momentum paired FFLO state [H, [ii, [13, [IH . The de- 
tailed structure of the proposed FFLO state is not clear, 
even in the homogeneous situation. How to extend the 
current BdG formalism to incorporate the FFLO state is 
therefore a fascinating issue. Another possibility is that 
the extremely narrow window for FFLO states in 3D is 
closed or reduced in size due to the presence of the har- 
monic trap. 

We have reported various mean-field thermodynam- 
ical properties of the imbalanced Fermi gas at unitar- 
ity, which is believed to be qualitatively reliable at low 
temperature. The experimentally observed bimodal dis- 
tribution in the profile of the minority spin state has 
been reproduced. We have determined the BCS super- 
fluid transition temperature as a function of the popu- 
lation imbalance, and have shown that it is consistent 
with recent experimental measurements. Quantitative 
theories of the imbalanced Fermi gas that take account 
of large quantum fluctuations that occur in the strongly 
correlated unitarity regime still need to be developed. A 
promising approach is to take into account pair fluctua- 
tions within the ladder approximation jg^l • This problem 
will be addressed in a future publication. 



VI. CONCLUSIONS 

We have developed an efficient and accurate hybrid 
procedure to solve the mean-fleld BdG equations for a 
spherically trapped Fermi gas with spin population im- 
balance. This enables us to thoroughly examine the ex- 
tensively used LDA approach. For a moderate large par- 
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